#ifndef BASEUTILS_H
#define BASEUTILS_H

#include <fstream>

#include "Math/Functor.h"
#include "Math/Factory.h"

#include "TASImage.h"
#include "TCanvas.h"
#include "TColor.h"
#include "TChain.h"
#include "TF1.h"
#include "TFile.h"
#include "TGaxis.h"
#include "TGraph.h"
#include "TGraphAsymmErrors.h"
#include "TGraphErrors.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TLine.h"
#include "TMath.h"
#include "TMatrixD.h"
#include "TMinuit.h"
#include "TLatex.h"
#include "TLegend.h"
#include "TPaletteAxis.h"
#include "TRandom3.h"
#include "TROOT.h"
#include "TStopwatch.h"
#include "TStyle.h"
#include "TTree.h"
#include "TVector3.h"

#include "TTreeStream.h"

#define EPSILON 1e-12


#ifndef TMATRIXDSYM_H
#include "TMatrixDSym.h"
#endif

class TAxis;
class TCanvas;
class TChain;
class TF1;
class TGraph;
class TGraphAsymmErrors;
class TGraphErrors;
class TH1;
class TH1D;
class TH2D;
class TLatex;
class TLegend;
class TMinuit;
class TPad;
class TTree;

class BaseUtils
{
 public:
  
 private:
  BaseUtils(const BaseUtils &p);
  BaseUtils & operator=(const BaseUtils &p);

  //===================================================================
  //                                Functions
  //===================================================================
 public:
  typedef Double_t (*FFunc)(const Double_t *xx, const Double_t *par);

  static Double_t ProtectedExp(const Double_t xx);
  static Double_t ProtectedPow(const Double_t xx, const Double_t yy);

  static Double_t ResolutionGaus(const Double_t *xx, const Double_t *par);

  static Double_t GausN(const Double_t *xx, const Double_t *par, const Int_t ntype, FFunc basefunc, const Int_t npartype);

  //===================================================================
  //                               Algorithm
  //===================================================================
 public:
  static Double_t GetCorrelatedError(const Double_t v0, const Double_t v1, const Double_t a0, const Double_t a1, const Double_t a2);

  static void Polint(const Int_t n, const Double_t *xa, const Double_t *ya, const Double_t x, Double_t &y, Double_t &dy, const Int_t kprint=0);

  //===================================================================
  //                                Fitting
  //===================================================================
 public:

  enum opt{
    kTMINUIT=0
  };

  static TMatrixDSym * GetCov(TMinuit * tmn, const Bool_t ifixed[], const Double_t errs[]);

  static TH1D * GetHfit(const TString hname, FFunc func, const Double_t par[], const Double_t xmin, const Double_t xmax, const Bool_t kbinlogx=0);

  static Int_t BinnedLikelihoodFit(const TH1 *hh, FFunc ffunc, const Int_t npar, Double_t pars[], Double_t errs[]=0x0, const Bool_t *pfix=0x0, Double_t *cov=0x0, const Double_t *lows=0x0, const Double_t *highs=0x0);

  static Int_t ChisquareFit(const Int_t nx, const Double_t xdata[], const Double_t ydata[], const Double_t *xerr, const Double_t yerr[], FFunc ffunc, const Int_t npar, Double_t pars[], Double_t *errs=0x0, Double_t *chi=0x0, const Bool_t *pfix=0x0, Double_t *cov=0x0, const Double_t *lows=0x0, const Double_t *highs=0x0);

  static void SetFitPrintLevel(const Int_t lev){ fgFitPrintLevel = lev; }
  static void SetTol(const Double_t tol){fTol = tol;}
  static void SetMaxIter(const Int_t iter){fMaxIter = iter;}
  static void SetMinOpt(const Int_t opt){fMinOpt=opt;}

 private:
  static Int_t FitKernel(const Int_t kml,const Int_t npar, Double_t pars[], Double_t errs[], Double_t chi[], const Bool_t *pfix, Double_t *cov, const Double_t *lows, const Double_t *highs);
  static Int_t TMinuitFitKernel(const Int_t kml,const Int_t npar, Double_t pars[], Double_t errs[], Double_t chi[], const Bool_t *pfix, Double_t *cov, const Double_t *lows, const Double_t *highs);

  static Int_t BinnedFit(const Bool_t kml, const Int_t nx, const Double_t xdata[], const Int_t nbin, FFunc ffunc, const Int_t npar, Double_t pars[], Double_t errs[], Double_t chi[], const Double_t thres, const Bool_t kxerr, const Bool_t *pfix, Double_t *cov, const Double_t *lows, const Double_t *highs);
  static Int_t BinnedFit(const Bool_t kml, const TH1 *hh, FFunc ffunc, const Int_t npar, Double_t pars[], Double_t errs[], Double_t chi[], const Double_t thres, const Bool_t kxerr, const Bool_t *pfix, Double_t *cov, const Double_t *lows, const Double_t *highs);
  static Int_t BinnedFit(const Bool_t kml, const Int_t nx, const Double_t xdata[], const Double_t ydata[], const Double_t *xerr, const Double_t yerr[], FFunc ffunc, const Int_t npar, Double_t pars[], Double_t errs[], Double_t chi[], const Bool_t *pfix,Double_t *cov, const Double_t *lows, const Double_t *highs);

  static Double_t GetFCNChisquare(const Double_t pars[]);
  static void MinFCNChisquare(Int_t &npar, Double_t *gin, Double_t &f, Double_t *pars, Int_t iflag);
  static Double_t GetFCNBinnedLikelihood(const Double_t pars[]);
  static void MinFCNBinnedLikelihood(Int_t &npar, Double_t *gin, Double_t &f, Double_t *pars, Int_t iflag);

  static void BinLogX(TAxis *axis);
  //=================

  static Int_t fgFitPrintLevel;  //print level

  static Int_t fgFitNStep; //n step
  static Int_t fgFitNData; //n data
  static Double_t *fgFitX; //data x
  static Double_t *fgFitY; //data y
  static Double_t *fgFitEx; //err x
  static Double_t *fgFitEy; //err y
  static FFunc fgFitFunc; //fit func
  static Double_t fTol;
  static Int_t fMaxIter;
  static Int_t fMinOpt;
  //===================================================================
  //                               Histograms
  //===================================================================
 public:
  static TH1D* ToPDF(const TH1 *h0, const TString hn);
  static TH2D* NormalHist(const TH2D *h0, const Double_t thres=0, const Bool_t kmax=kFALSE);

  static Double_t * GetAxisArray(TAxis * aa);

  static void FitSlicesY(const TH2D *hh, FFunc ffunc, TH1D *& hnor, TH1D *&hmpv, TH1D *&hwid, TH1D *&hres, TH1D *&hchi, const Double_t thres=10, const Double_t lowfrac=0.0, const Double_t highfrac=1.0);

  static TGraphAsymmErrors * TreeToGraph(const TString tname, const TString cut, Char_t * sx, Char_t * sy, Char_t * sex0=(Char_t*)"", Char_t * sex1=(Char_t*)"", Char_t * sey0=(Char_t*)"", Char_t * sey1=(Char_t*)"", const Bool_t kXLogToLinear=kFALSE, const Bool_t kprint=kFALSE);

  //===================================================================
  //                           IO
  //===================================================================
 public:
  static TChain * InputFiles(const TString file, const TString tr, Char_t * dir=0x0);

  //===================================================================
  //
  //===================================================================
};

Int_t BaseUtils::fgFitPrintLevel = 1;
Int_t BaseUtils::fgFitNStep=0;
Int_t BaseUtils::fgFitNData=0;
Double_t * BaseUtils::fgFitX = 0x0;
Double_t * BaseUtils::fgFitY = 0x0;
Double_t * BaseUtils::fgFitEx = 0x0;
Double_t * BaseUtils::fgFitEy = 0x0;
BaseUtils::FFunc BaseUtils::fgFitFunc = 0x0;
Double_t BaseUtils::fTol = 1;
Int_t BaseUtils::fMaxIter = 1e8;
Int_t BaseUtils::fMinOpt = BaseUtils::kTMINUIT;

//===================================================================  
//                               Algorithm 
//===================================================================
Double_t BaseUtils::GetCorrelatedError(const Double_t v0, const Double_t v1, const Double_t a0, const Double_t a1, const Double_t a2)
{
  //                   ( a0   a1 ) (v0)
  //df^2 = (v0,v1) *   ( a1   a2 ) (v1) = v0*v0*a0 + v1*v1*a2 + 2*v0*v1*a1
  //f=f(a,b), v0=df/da, v1=df/db
  //f=a/b, v0=f/a, v1=-f/b

  const Double_t df2 = v0*v0*a0 + v1*v1*a2 + 2*v0*v1*a1;
  /*
  //BaseUtils::GetCorrelatedError df2<0 -2.980232e-08 1.000000e-10 -- 1.640908e+07 -1.640908e+07 -- 3.221823e-07 3.221823e-07 3.221823e-07
  const Double_t eps = 1e-2;
  if(df2<eps){
    if(TMath::Abs(df2)<eps){
      return 0;
    }
    else{
      printf("BaseUtils::GetCorrelatedError df2<0 %e %e -- %e %e -- %e %e %e\n", df2, eps, v0, v1, a0, a1, a2);exit(1);
    }
  }
  */
  if(df2<EPSILON){
    return 0;
  }
  return TMath::Sqrt( df2 );
}

//Numerical Recipes in C, 1997, 3.1 Polynomial Interpolation and Extrapolation, p109
//Neville's algorithm
//error defined here depends on the choice of the route from ya[ns] to y_final
void BaseUtils::Polint(const Int_t n, const Double_t *xa, const Double_t *ya, const Double_t x, Double_t &y, Double_t &dy, const Int_t kprint)
{
  /*
  for(Int_t ii=0; ii<n; ii++){
    printf("BaseUtils::Polint %10.10e %03d/%03d: %10.10e %10.10e\n", x, ii, n, xa[ii], ya[ii]);
  }
  */

  Double_t c[n];
  Double_t d[n];
  Double_t dift=-999;

  Int_t ns=0;
  Double_t dif = fabs(x-xa[ns]); 
  for (Int_t i=0; i<n;i++) {
    if ( (dift=fabs(x-xa[i])) < dif) { 
      ns=i; 
      dif=dift; 
    }
    c[i]=ya[i]; 
    d[i]=ya[i];

    if(kprint) printf("polint =====> x[%d]: %3.10e\ty[%d]: %3.10e\n", i, xa[i], i, ya[i]);
  }

  y = ya[ns];
  if(kprint) printf("======> y0: %3.10e at ns=%d\n", y, ns);

  Int_t did=ns-1;
  for (Int_t m=1;m<=n-1;m++){
    for (Int_t i=0; i<n-m; i++){
      Double_t ho=xa[i]-x; 
      Double_t hp=xa[i+m] -x;
      Double_t w = c[i+1]-d[i];
      Double_t den=ho-hp;
      if ( den == 0.0){
        printf("BaseUtils::Polint error (den=ho-hp) == 0.0!! %e %e %e\n", x, xa[i], xa[i+m]); exit(1);
      }

      c[i]=ho*w/den;
      d[i]=hp*w/den; 
      //if(kprint) printf("======> m %d, c[%d]: %5.6e\t d[%d]: %5.6e\n", m,i,c[i],i, d[i]);
    }
    if( m <= n -1-ns ){
      dy = c[ns];
      if(kprint) printf("------------> m %d, c[%d]: %5.6e\n", m, ns, dy);
    }
    else {
      dy = d[did];
      if(kprint) printf("------------> m %d, d[%d]: %5.6e\n", m, did, dy);
      did--;
    }
    y += dy;
    if(kprint) printf("======> y%d: %3.10e\n",m, y);
  }
}

//===================================================================
//                                Functions
//===================================================================
Double_t BaseUtils::ProtectedExp(const Double_t xx)
{
  /*
    root [14] exp(709)
    (const double)8.21840746155497238e+307
    root [15] exp(710)
    Error: exp param[0]=710 up:709 low:-inf out of range (tmpfile):1:
  */
  const Double_t maxexparg = 700;
  
  return TMath::Exp(TMath::Min(maxexparg, xx));
}

Double_t BaseUtils::ProtectedPow(const Double_t xx, const Double_t yy)
{
  if(xx<0){
    printf("BaseUtils::ProtectedPow xx<0 %e\n", xx); exit(1);
  }
  //xx^yy = exp(log(xx^yy))=exp(yy log(xx))
  return ProtectedExp(yy*TMath::Log(xx+EPSILON));
}

Double_t BaseUtils::ResolutionGaus(const Double_t *xx, const Double_t *par)
{
  //npar = 3

  //
  //root [0] TMath::Gaus(3, 0, -9,1)
  //(Double_t)(-4.19314697436659056e-02)
  //root [1] TMath::Gaus(3, 0, 9,1)
  //(Double_t)4.19314697436659056e-02
  //
  //!! width has to be abs-ed!!
  return par[0]*TMath::Gaus(xx[0], par[1], par[1]*TMath::Abs(par[2]), kTRUE);
}

Double_t BaseUtils::GausN(const Double_t *xx, const Double_t *par, const Int_t ntype, FFunc basefunc, const Int_t npartype)
{
  //npar = ntype*npartype

  Double_t sum = 0;

  for(Int_t itype=0; itype<ntype; itype++){
    sum += basefunc(xx, &(par[npartype*itype]));
  }

  return sum;
}

//===================================================================
//                                Fitting 
//===================================================================
TMatrixDSym * BaseUtils::GetCov(TMinuit * tmn, const Bool_t ifixed[], const Double_t errs[])
{
  // case of fixed params need to take care
  //http://root.cern.ch/root/html/src/TMinuitMinimizer.cxx.html#Ilx.Z

  const Int_t nfree = tmn->GetNumFreePars();
  const Int_t npar = tmn->GetNumPars();

  TMatrixDSym matfree(nfree); 
  tmn->mnemat(matfree.GetMatrixArray(), nfree);

  TMatrixDSym *matall = new TMatrixDSym(npar);

  //r: row, c: column
  Int_t rfree = 0;
  for(Int_t rr=0; rr<npar; rr++){
    if(ifixed[rr]){
      continue;
    }
    Int_t cfree = 0;
    for(Int_t cc=0; cc<npar; cc++){
      if(ifixed[cc]){
        continue;
      }
      (*matall)[rr][cc] = matfree[rfree][cfree];
      cfree++;
    }
    rfree++;
  }

  /*
  //test->
    for(int ii=0; ii<NparBin(); ii++){
    const Int_t id = NparBuffer()+ii;
    printf("ii %d: err %e %e -- %e\n", ii, fLowErrs[id], sqrt((*matall)[id][id]), fLowErrs[id]-sqrt((*matall)[id][id]));
    for(int jj=0; jj<NparBin(); jj++){
    const Int_t i2 = NparBuffer()+jj;
    printf("%e ", (*matall)[id][i2]);
    }
    printf("\n");
    }
    //<-test
  */

  //check consistency!
  const Double_t eps = 1e-1;
  for(Int_t ipar=0; ipar<npar; ipar++){
    const Double_t cov = (*matall)[ipar][ipar];
    const Double_t sigma = TMath::Sqrt(cov);
    if(TMath::Abs(errs[ipar]-sigma)/errs[ipar]>eps){
      printf("BaseUtils::GetCov warning!! --> err-sigma > eps %d: %e %e %e -- %e\n", ipar, errs[ipar], cov, sigma, eps);
    }
  }

  return matall;
}

Int_t BaseUtils::FitKernel(const Int_t kml,const Int_t npar, Double_t pars[], Double_t errs[], Double_t chi[], const Bool_t *pfix, Double_t *cov, const Double_t *lows, const Double_t *highs)
{
  //
  //kml: 0: chi2, 1: maximum likelihood, 2: binned likelihood
  //
  static Int_t counter = 0;
  if(counter==0){
    printf("\n\n**************** BaseUtils::FitKernel fMinOpt %d **************** \n", fMinOpt);
  }

  counter++;
  if(fMinOpt!=kTMINUIT){
    printf("non TMinuit not implemented!\n");exit(1);
  }

  return TMinuitFitKernel(kml, npar, pars, errs, chi, pfix, cov, lows, highs);
}


Int_t BaseUtils::TMinuitFitKernel(const Int_t kml,const Int_t npar, Double_t pars[], Double_t errs[], Double_t chi[], const Bool_t *pfix, Double_t *cov, const Double_t *lows, const Double_t *highs)
{
  //
  //fitkernel for both ML and chi2 fit
  //

  TStopwatch stwa;
  if(fgFitPrintLevel>=0){
    stwa.Start(kTRUE);
  }
  //====>

  TMinuit * mnt = new TMinuit(npar);
  switch(kml){
  case 0:
    mnt->SetFCN(MinFCNChisquare);
    break;
  case 2:
    mnt->SetFCN(MinFCNBinnedLikelihood);
    break;
  default:
    printf("BaseUtils::TMinuitFitKernel bad kml %d\n", kml);exit(1);
  }
  mnt->SetPrintLevel(fgFitPrintLevel);
    
  //http://wwwasdoc.web.cern.ch/wwwasdoc/WWW/minuit/minmain/chapter2_8.html#SECTION0081000000000000000
  //SET ERRordef<up>
  //Sets the value of UP (default value= 1.), defining parameter errors. Minuit defines parameter errors as the change in parameter value required to change the function value by UP. Normally, for chisquared fits UP=1, and for negative log likelihood, UP=0.5. 
  const Double_t errdef = kml? 0.5 : 1;
  mnt->SetErrorDef(errdef);
    
  Int_t nfree = 0;
  for(Int_t ipar=0;ipar<npar;ipar++){
    const Double_t start=pars[ipar];
    //1e-3*par seems too small
    const Double_t ss= TMath::Abs(pars[ipar])>EPSILON ? TMath::Abs(pars[ipar])*6e-1 : 1e-3;//important if limits are set
    //1e-2; //restore to old values for TRD performance /u/xlu/.task/dedx/drawSummary in /u/xlu/.task/dedx/AllDataSummary/hist/draw_UpdatedStyle_Comb
    //
    const Double_t ll=(lows && lows[ipar]>-1e9)?lows[ipar]:0;
    const Double_t hh=(highs && highs[ipar]>-1e9)?highs[ipar]:0;
    
    mnt->DefineParameter(ipar, Form("p%d",ipar), start, ss, ll,hh);
    nfree++;
    if(pfix && pfix[ipar]){
      mnt->FixParameter(ipar);
      nfree--;
    }
  }

  //http://wwwasdoc.web.cern.ch/wwwasdoc/WWW/minuit/minmain/chapter2_8.html#SECTION0081000000000000000
  //MIGrad[maxcalls] [tolerance]
  //Causes minimization of the function by the method of Migrad, the most efficient and complete single method, recommended for general functions (see also MINImize). The minimization produces as a by-product the error matrix of the parameters, which is usually reliable unless warning messages are produced. The optional argument [maxcalls] specifies the (approximate) maximum number of function calls after which the calculation will be stopped even if it has not yet converged. The optional argument [tolerance] specifies required tolerance on the function value at the minimum. The default tolerance is 0.1, and the minimization will stop when the estimated vertical distance to the minimum (EDM) is less than 0.001*[tolerance]*UP (see [SET ERRordef]SET ERR). 

  Double_t arg[2]={fMaxIter, fTol};
  Int_t migradflag = -999;
  mnt->mnexcm("MIGRAD", arg, 2, migradflag);

  //get par right after MIGRAD before HESSE modifies them, initial pars not changed
  Double_t migp[10];
  Double_t dummy=-999;
  for(Int_t ipar=0;ipar<npar;ipar++){
    mnt->GetParameter(ipar,migp[ipar], dummy);
  }

  Int_t hesseflag = -999;
  mnt->mnexcm("HESSE", arg ,1, hesseflag);

  const Int_t kfail = migradflag+hesseflag;
  if(kfail){
    printf("BaseUtils::FitKernel fitting error !! %d %d keep raw par and set err to -999!!!\n", migradflag, hesseflag);
    if(errs){
      for(Int_t ipar=0;ipar<npar;ipar++){
        errs[ipar]=-999;
      }
    }
    if(chi){
      chi[0]=chi[1]=-999;
    }
  }
  else{
    for(Int_t ipar=0;ipar<npar;ipar++){
      if(errs) {
        mnt->GetParameter(ipar, dummy,errs[ipar]);
      }

      pars[ipar]=migp[ipar];
    }
    if(chi){
      chi[0] = GetFCNChisquare(pars);
      if(nfree!=mnt->GetNumFreePars()){
        printf("BaseUtils::TMinuitFitKernel nfree!=mnt->GetNumFreePars() %d %d\n", nfree, mnt->GetNumFreePars());exit(1);
      }
      chi[1] = -nfree;//need to add nx outside FitKernel
    }
  }

  if(cov){
    mnt->mnemat(cov,npar);
  }

  delete mnt;

  //<===
 
  if(fgFitPrintLevel>=0){
    stwa.Stop(); 
    stwa.Print("m"); //"m": milli sec; "u": micro sec
  } 
 
  //test->
  /*
  for(Int_t ii=0; ii<npar; ii++){
    printf("test0 par %d %e\n", ii, pars[ii]);
    printf("test1 err %d %e\n", ii, errs[ii]);
  }
  */
  //<-test

  return kfail;
}


Int_t BaseUtils::BinnedFit(const Bool_t kml, const Int_t nx, const Double_t xdata[], const Int_t nbin, FFunc ffunc, const Int_t npar, Double_t pars[], Double_t errs[], Double_t chi[], const Double_t thres, const Bool_t kxerr, const Bool_t *pfix, Double_t *cov, const Double_t *lows, const Double_t *highs)
{
  //
  //chi2fit from unbinned X -> TH1 
  //
  Int_t index[nx];
  TMath::Sort(nx, xdata, index, kFALSE);
  const Double_t xmin = xdata[index[0]];
  const Double_t xmax = xdata[index[nx-1]];
  const Double_t dx = (xmax-xmin)/nbin;

  TH1D hh("chisquarefittmphh","",nbin, xmin-dx, xmax+dx);
  for(Int_t ii=0; ii<nx; ii++){
    hh.Fill(xdata[ii]);
  }

  while(hh.GetBinContent(hh.GetMaximumBin())<10 && hh.GetNbinsX()>10){
    hh.Rebin(2);
  }

  Double_t tmperr[npar];
  Double_t tmpch[10];
  const Int_t kfail = BinnedFit(kml, &hh, ffunc, npar, pars, tmperr, tmpch, thres, kxerr, pfix, cov,lows, highs);
  if(errs){
    memcpy(errs, tmperr, sizeof(tmperr));
  }
  if(chi){
    memcpy(chi, tmpch, sizeof(Double_t)*2);
  }
  return kfail;
}


Int_t BaseUtils::BinnedLikelihoodFit(const TH1 *hh, FFunc ffunc, const Int_t npar, Double_t pars[], Double_t errs[], const Bool_t *pfix, Double_t *cov, const Double_t *lows, const Double_t *highs)
{
  //keep empty bins
  const Double_t thres = -1;
  return BinnedFit(kTRUE, hh, ffunc, npar, pars, errs, 0x0, thres, kFALSE, pfix, cov, lows, highs);
}

Int_t BaseUtils::BinnedFit(const Bool_t kml, const TH1 *hh, FFunc ffunc, const Int_t npar, Double_t pars[], Double_t errs[], Double_t chi[], const Double_t thres, const Bool_t kxerr, const Bool_t *pfix, Double_t *cov, const Double_t *lows, const Double_t *highs)
{
  //
  //chisquare fit from TH1
  //

  const Int_t nmax = hh->GetNbinsX();
  Double_t xdata[nmax];
  Double_t ydata[nmax];
  Double_t xerr[nmax];
  Double_t yerr[nmax];

  Int_t nx=0;
  for(Int_t ii=hh->GetXaxis()->GetFirst(); ii<=hh->GetXaxis()->GetLast(); ii++){
    const Double_t ey = hh->GetBinError(ii);
    const Double_t yy = hh->GetBinContent(ii);
    //keep all for kml==kTRUE
    //very important for BinnedLikelihood!! only keeping empty bins can the result be good
    if(!kml){
      //skip empty bin
      if(ey<EPSILON){
        if(yy>EPSILON){
          printf("BaseUtils::ChisquareFit error! %d %e %e\n", ii, ey, yy); exit(1);
        }
        continue;
      }

      //skip low statstics bin
      if(yy<thres){
        continue;
      }
    }
    const Double_t xx = hh->GetBinCenter(ii);
    const Double_t ex = hh->GetBinWidth(ii)/TMath::Sqrt(12);

    xdata[nx]=xx;
    ydata[nx]=yy;
    xerr[nx]=ex;
    yerr[nx]=ey;
    nx++;
  }

  return BinnedFit(kml, nx, xdata, ydata, kxerr?xerr:0x0, yerr, ffunc, npar, pars, errs, chi, pfix, cov, lows, highs);
}

Int_t BaseUtils::ChisquareFit(const Int_t nx, const Double_t xdata[], const Double_t ydata[], const Double_t *xerr, const Double_t yerr[], FFunc ffunc, const Int_t npar, Double_t pars[], Double_t errs[], Double_t chi[], const Bool_t *pfix, Double_t *cov, const Double_t *lows, const Double_t *highs)
{
  return BinnedFit(kFALSE, nx, xdata, ydata, xerr, yerr, ffunc, npar, pars, errs, chi, pfix, cov, lows, highs);
}

Int_t BaseUtils::BinnedFit(const Bool_t kml, const Int_t nx, const Double_t xdata[], const Double_t ydata[], const Double_t *xerr, const Double_t yerr[], FFunc ffunc, const Int_t npar, Double_t pars[], Double_t errs[], Double_t chi[], const Bool_t *pfix, Double_t *cov, const Double_t *lows, const Double_t *highs)
{
  //
  //routine for chisquare fit from arrays X and Y
  //
  //
  if(nx<=0){
    printf("BaseUtils::ChisquareFit error nx<=0 %d\n", nx); exit(1);
  }

  fgFitNStep=0;
  fgFitNData = nx;
  fgFitX = new Double_t[nx];
  fgFitY = new Double_t[nx];
  fgFitEx = new Double_t[nx];
  fgFitEy = new Double_t[nx];
  for(Int_t ix=0; ix<fgFitNData; ix++){
    fgFitX[ix]=xdata[ix];
    fgFitY[ix]=ydata[ix];
    fgFitEx[ix]=xerr?xerr[ix]:0;
    fgFitEy[ix]=yerr[ix];
    if(fgFitPrintLevel>=0){
      printf("BaseUtils::ChisquareFit %d: %e %e %e %e\n", ix, fgFitX[ix], fgFitY[ix], fgFitEx[ix], fgFitEy[ix]);
    }
  }
  fgFitFunc = ffunc;

  const Int_t kfail = FitKernel(kml?2:0, npar, pars, errs, chi, pfix, cov, lows, highs);

  if(chi){
    chi[1] += nx;
    printf("\tBaseUtils::ChisquareFit Chi2/NDOF %.1f/%.0f = %.1f\n\n", chi[0], chi[1], chi[1]? chi[0]/chi[1] : -999);
  }

  delete fgFitX;
  delete fgFitY;
  delete fgFitEx;
  delete fgFitEy;

  return kfail;
}

Double_t BaseUtils::GetFCNChisquare(const Double_t pars[])
{
  //
  //chi2
  //

  //const Int_t ktest = 1000000;//207;//4165;

  Double_t f = 0;

  /*
  if(fgFitNStep>=ktest){
    for(Int_t ii=0; ii<npar; ii++){
      printf("in MinFCN pars %d %e\n", ii, pars[ii]);
    }
  }
  */

  for(Int_t idata=0; idata<fgFitNData; idata++){
    //if(fgFitNStep>=ktest) printf("test %d %e %e %e %e %e %e\n", idata, fgFitX[idata], pars[0], pars[1], pars[2], pars[3], pars[4]);

    const Double_t evaly = fgFitFunc(&(fgFitX[idata]), pars);
    const Double_t dy = fgFitY[idata] - evaly;
    
    Double_t dfdx=0;
    //calculated df/dx if necessary
    if( fgFitEx[idata]!=0){
      const Double_t derix0=fgFitX[idata]-fgFitEx[idata];
      const Double_t derix1=fgFitX[idata]+fgFitEx[idata];
      const Double_t deriy0=fgFitFunc(&derix0, pars);
      const Double_t deriy1=fgFitFunc(&derix1, pars);
      dfdx = (deriy1-deriy0)/(derix1-derix0);
    }
    f += dy*dy/(fgFitEy[idata]*fgFitEy[idata] + dfdx*fgFitEx[idata]*dfdx*fgFitEx[idata]);
    
    //if(fgFitNStep>=ktest) printf("in MinFCN idata %d p0 %e p1 %e data %e f %e dy %e ey %e\n", idata, pars[0], pars[1], fgFitX[idata], f, dy, fgFitEy[idata]);
  }

  //if(fgFitNStep>=ktest) printf("fgFitNStep %d f %e\n", fgFitNStep, f);
  fgFitNStep++;

  return f;
}

void BaseUtils::MinFCNChisquare(Int_t &npar, Double_t *gin, Double_t &f, Double_t *pars, Int_t iflag)
{
  //
  //function to be minimized in chi2 fit
  //

  f = GetFCNChisquare(pars);
}

Double_t BaseUtils::GetFCNBinnedLikelihood(const Double_t pars[])
{
  //
  //binned likelihood
  //

  //const Int_t ktest = 0;//1000000;//207;//4165;

  //the normalization is important for the error estimation, but not the mean!
  //the norm should be the same for all bins, so just get it from first non-0 bin
  Double_t norm = 0;
  for(Int_t idata=0; idata<fgFitNData; idata++){
    if(fgFitEy[idata]>EPSILON){
      if(fgFitY[idata]<EPSILON){
        printf("BaseUtils::GetFCNBinnedLikelihood error y<eps, while ey>eps!! %e %e %e", EPSILON,fgFitEy[idata],fgFitY[idata] );exit(1);
      }
      //y = kN, ey = k sqrt(N), y/ey = sqrt(N)
      const Double_t trueCount = TMath::Power(fgFitY[idata]/fgFitEy[idata],2);
      norm = trueCount/fgFitY[idata];
      break;
    }
  }
  if(norm<EPSILON){
    printf("BaseUtils::GetFCNBinnedLikelihood error norm 0 !! %e\n", norm); exit(1);
  }

  Double_t f = 0;

  for(Int_t idata=0; idata<fgFitNData; idata++){
    const Double_t evaly = fgFitFunc(&(fgFitX[idata]), pars);
    const Double_t prob = TMath::Poisson(fgFitY[idata]*norm, evaly*norm);
    f += -TMath::Log( prob>EPSILON ? prob : EPSILON);
    
    //if(fgFitNStep>=ktest) printf("in MinFCN idata %d p0 %e p1 %e p2 %e x %e y %e y0 %e prob %e f %e\n", idata, pars[0], pars[1], pars[2], fgFitX[idata], fgFitY[idata], evaly, prob, f);
  }

  //if(fgFitNStep>=ktest) printf("fgFitNStep %d f %e\n", fgFitNStep, f);
  //if(fgFitNStep>=ktest) printf("test\n"); exit(1);

  fgFitNStep++;
  return f;
}

void BaseUtils::MinFCNBinnedLikelihood(Int_t &npar, Double_t *gin, Double_t &f, Double_t *pars, Int_t iflag)
{
  //
  //function to be minimized in chi2 fit
  //

  f = GetFCNBinnedLikelihood(pars);
}


void BaseUtils::BinLogX(TAxis *axis)
{
  //
  // Method for the correct logarithmic binning of histograms
  // copied and modified from AliTPCcalibBase

  const Int_t bins = axis->GetNbins();

  const Double_t from = axis->GetXmin();
  const Double_t to = axis->GetXmax();
  if (from<EPSILON){
    //printf("BaseUtils::BinLogX warning xmin < epsilon! nothing done, axis not set. %e\n", from);  exit(1);                                                                                                                                       
    return;
  }

  Double_t *new_bins = new Double_t[bins + 1];

  new_bins[0] = from;
  const Double_t factor = pow(to/from, 1./bins);

  for (int i = 1; i <= bins; i++) {
    new_bins[i] = factor * new_bins[i-1];
  }
  axis->Set(bins, new_bins);
  delete [] new_bins;
}


TH1D * BaseUtils::GetHfit(const TString hname, FFunc func, const Double_t par[], const Double_t xmin, const Double_t xmax, const Bool_t kbinlogx)
{
  TH1D * h1=new TH1D(hname, "", 500, xmin, xmax);
  if(kbinlogx){
    BinLogX(h1->GetXaxis());
  }
  for(Int_t ii=1; ii<=h1->GetNbinsX(); ii++){
    const Double_t xx = h1->GetBinCenter(ii);
    const Double_t yy = func(&xx, par);
    h1->SetBinContent(ii, yy);
  }
  return h1;
}

//===================================================================
//                           Histograms
//===================================================================

TGraphAsymmErrors * BaseUtils::TreeToGraph(const TString tname, const TString cut, Char_t * sx, Char_t * sy, Char_t * sex0, Char_t * sex1, Char_t * sey0, Char_t * sey1, const Bool_t kXLogToLinear, const Bool_t kprint)
{
  TTree * rawtree=(TTree*)gDirectory->Get(tname);
  if(!rawtree){
    printf("BaseUtils::TreeToGraph no tree %s\n", tname.Data());exit(1);
  }

  TTree *cuttree = rawtree->CopyTree(cut);
  TTree * tmptree=(TTree*) cuttree->Clone(Form("tmptree%s%s",tname.Data(), sy));//necessary to avoid address confusion when executed repeatedly.

  Double_t vx = -1e10, vy = -1e10, vex0 =-1e10, vex1 = -1e10, vey0 = -1e10, vey1 = -1e10;
  tmptree->SetBranchAddress(sx, &vx);
  tmptree->SetBranchAddress(sy, &vy);
  if(strlen(sex0) && strlen(sex1)){
    tmptree->SetBranchAddress(sex0, &vex0);
    tmptree->SetBranchAddress(sex1, &vex1);
  }
  if(strlen(sey0)){
    tmptree->SetBranchAddress(sey0, &vey0);
  }
  //if sey0 and sey1 are the same, vey0 will not be set from tree, has to be set by vey1
  if(strlen(sey1)){
    tmptree->SetBranchAddress(sey1, &vey1);
  }
  
  TGraphAsymmErrors * tmpgr = new TGraphAsymmErrors;
  const Int_t ntot = tmptree->GetEntries();
  for(Int_t ii=0; ii<ntot; ii++){
    tmptree->GetEntry(ii);

    //remove trivial -1e10
    if(vex0<-1e9) vex0=0;
    if(vex1<-1e9) vex1=0;
    if(vey0<-1e9) vey0=0;
    if(vey1<-1e9) vey1=0;

    //minos error low < 0 by convention
    vex0=TMath::Abs(vex0);
    vex1=TMath::Abs(vex1);
    vey0=TMath::Abs(vey0);
    vey1=TMath::Abs(vey1);

    if(strcmp(sey0,sey1)==0){
      vey0=vey1;
    }

    if(kXLogToLinear){
      const Double_t xtrue = pow(10, vx);
      const Double_t xlow = pow(10, vx-vex0);
      const Double_t xhigh = pow(10, vx+vex1);
      vx = xtrue;
      vex0 = xtrue - xlow;
      vex1 = xhigh - xtrue;
    }

    tmpgr->SetPoint(ii, vx, vy);
    tmpgr->SetPointError(ii, vex0, vex1, vey0, vey1);
    if(kprint){
      printf("BaseUtils::TreeToGraph %s %s %s -- %d %e %e -- %e %e %e %e\n", tname.Data(), sx, sy, ii, vx, vy,  vex0, vex1, vey0, vey1);
    }
  }
  delete tmptree;

  TGraphAsymmErrors * gr = (TGraphAsymmErrors *)tmpgr->Clone(Form("gr%s%s",tname.Data(), sy));
  delete tmpgr;
  return gr;
}

TH1D* BaseUtils::ToPDF(const TH1 *hraw, const TString hn)
{
  TH1D * hist = (TH1D*) hraw->Clone((hn+hraw->GetName())+"pdf");
  Double_t ntot = 0;
  hist->Scale(0);

  //only consider region in range, ignore over/underflow. Important for fitting experiment data where under/overflow is ignored.
  const Int_t x0 = hist->GetXaxis()->GetFirst();
  const Int_t x1 = hist->GetXaxis()->GetLast();

  for(Int_t ii=x0; ii<=x1; ii++){
    const Double_t err = hraw->GetBinError(ii);
    const Double_t cont = hraw->GetBinContent(ii);

    //skip empty bins
    if(err<EPSILON){
      if(cont>EPSILON){
        printf("BaseUtils::ToPDF error! %d %e %e\n", ii, err, cont); exit(1);
      }
      continue;
    }

    const Double_t nn = cont*cont/err/err;
    hist->SetBinContent(ii,nn);
    hist->SetBinError(ii, TMath::Sqrt(nn));
    ntot += nn;
  }

  if(ntot<EPSILON){
    printf("BaseUtils::ToPDF ntot<epsilon ! %f %f\n", hist->GetEntries(), ntot);
    exit(1);
  }
  
  const Double_t bw = hist->GetBinWidth(1);
  for(Int_t ib=x0; ib<=x1; ib++){
    const Double_t cont = hist->GetBinContent(ib);
    if(cont<EPSILON)
      continue;

    //in case of finit number of bins (i.e. eff not always small), Binomial error is more accurate than Poisson error
    const Double_t eff = cont/ntot;
    const Double_t pdf = eff/bw;

    const Double_t dpdf = sqrt(eff*(1-eff)/ntot) / bw;
    hist->SetBinContent(ib, pdf);
    hist->SetBinError(ib, dpdf);
  }

  hist->SetEntries(ntot);

  return hist;
}


TH2D* BaseUtils::NormalHist(const TH2D *hraw, const Double_t thres, const Bool_t kmax)
{
  //
  //normalized histogram
  //

  TH2D *hh=(TH2D*)hraw->Clone(Form("%snor",hraw->GetName()));
  hh->Scale(0);

  const Int_t x0 = hh->GetXaxis()->GetFirst();
  const Int_t x1 = hh->GetXaxis()->GetLast();
  const Int_t y0 = hh->GetYaxis()->GetFirst();
  const Int_t y1 = hh->GetYaxis()->GetLast();

  Double_t hmax = -1e10;
  Double_t hmin = 1e10;
  Double_t nent = 0;
  for(Int_t ix=x0; ix<=x1; ix++){

    //if option "e" is specified, the errors are computed. if option "o" original axis range of the taget axes will be kept, but only bins inside the selected range will be filled.
    TH1D * sliceh = hraw->ProjectionY(Form("tmpnormalhist%sx%d", hh->GetName(), ix), ix, ix, "oe");
    const Double_t tot = sliceh->GetEntries();

    TH1D * pdfh=0x0;

    if(tot>EPSILON){
      nent += tot;
    
      Double_t imax = -999;

      if(!kmax){
        pdfh = ToPDF(sliceh,"tmp");
      }
      else{
        imax = sliceh->GetBinContent(sliceh->GetMaximumBin());
      }

      for(Int_t iy=y0; iy<=y1; iy++){
        const Double_t cont = kmax ? sliceh->GetBinContent(iy)/imax : pdfh->GetBinContent(iy);
        const Double_t ierr = kmax ? sliceh->GetBinError(iy)/imax   : pdfh->GetBinError(iy);
        if(tot>thres && cont>0){
          hh->SetBinContent(ix, iy, cont);
          hh->SetBinError(ix,iy, ierr);
          if(cont>hmax) hmax = cont;
          if(cont<hmin) hmin = cont;
        }
      }
    }

    delete pdfh;
    delete sliceh;
  }

  hh->SetEntries(nent);
  hh->SetMinimum(0.99*hmin);
  hh->SetMaximum(1.1*hmax);
  return hh;
}

Double_t * BaseUtils::GetAxisArray(TAxis * aa)
{
  TArrayD * xs = new TArrayD(*(aa->GetXbins()));
  if(xs->GetSize()==0){
    const Int_t nbin=aa->GetNbins();
    xs->Set(nbin+1);
    for(Int_t ii=0; ii<=nbin; ii++){
      xs->SetAt(aa->GetBinUpEdge(ii),ii);
    }
  }

  Double_t * bins = xs->GetArray();
  if(!bins){
    printf("BaseUtils::GetAxisArray bins null!! %d\n", xs->GetSize());exit(1);
  }
  return bins;
}

void BaseUtils::FitSlicesY(const TH2D *hh, FFunc ffunc, TH1D *&hnor, TH1D *&hmpv, TH1D *&hwid, TH1D *&hres, TH1D *&hchi, const Double_t thres, const Double_t lowfrac, const Double_t highfrac)
{
  //
  //ffunc 3 parameter function: ALandau, AGaus
  //

  const Int_t x0 = hh->GetXaxis()->GetFirst();
  const Int_t x1 = hh->GetXaxis()->GetLast();
  const Int_t y0 = hh->GetYaxis()->GetFirst();
  const Int_t y1 = hh->GetYaxis()->GetLast();

  const Int_t nx = hh->GetNbinsX();
  const Int_t ny = hh->GetNbinsY();
  const Double_t xmin = hh->GetXaxis()->GetXmin();
  const Double_t xmax = hh->GetXaxis()->GetXmax();
  const Double_t ymin = hh->GetYaxis()->GetXmin();
  const Double_t ymax = hh->GetYaxis()->GetXmax();

  hnor = new TH1D(Form("%s_amp",hh->GetName()), "", nx, xmin, xmax);
  hmpv = new TH1D(Form("%s_mpv",hh->GetName()), "", nx, xmin, xmax);
  hwid = new TH1D(Form("%s_wid",hh->GetName()), "", nx, xmin, xmax);
  hres = new TH1D(Form("%s_res",hh->GetName()), "", nx, xmin, xmax);
  hchi = new TH1D(Form("%s_chi",hh->GetName()), "", nx, xmin, xmax);

  const Double_t *hxbins = GetAxisArray(hh->GetXaxis());

  hnor->GetXaxis()->Set(nx, hxbins);
  hmpv->GetXaxis()->Set(nx, hxbins);
  hwid->GetXaxis()->Set(nx, hxbins);
  hres->GetXaxis()->Set(nx, hxbins);
  hchi->GetXaxis()->Set(nx, hxbins);

  for(Int_t ix=x0; ix<=x1; ix++){
    TH1D *htmp = new TH1D(Form("FitSlicesY_%s_%d", hh->GetName(), ix),"",ny, ymin, ymax);
    //checked, ok
    const Double_t *hhybins = GetAxisArray(hh->GetYaxis());
    //hh->GetYaxis()->GetXbins()->GetArray();
    if(hhybins){
      htmp->GetXaxis()->Set(ny, hhybins);
    }

    Double_t ntot = 0;
    for(Int_t iy=y0; iy<=y1; iy++){
      const Double_t be = hh->GetBinError(ix,iy);
      const Double_t bc = hh->GetBinContent(ix, iy);

      if(be<EPSILON){
        if(bc>EPSILON){
          printf("BaseUtils::FitSlicesY error %d %d %e %e\n", ix, iy, be, bc); exit(1);
        }
        continue;
      }

      htmp->SetBinContent(iy, bc);
      htmp->SetBinError(iy, be);

      ntot += (bc/be)*(bc/be);

      //if(be) printf("test %d %d : %f %f %f\n", ix, iy, bc, be, pow(bc/be,2));
    }

    hnor->SetBinContent(ix, ntot);
    hnor->SetBinError(  ix, 0);
    
    if(ntot<thres || htmp->GetRMS()<EPSILON){
      delete htmp;
      continue;
    }

    //test htmp->Draw();
    Double_t pars[10]={htmp->Integral(0,htmp->GetNbinsX()+1)*htmp->GetBinWidth(1), htmp->GetMean(), htmp->GetRMS()};
    Double_t errs[10], chi[10];
    
    if(ffunc){
      printf("not implemented!\n");exit(1);
    }
    else{
      printf("truncated mean removed!!\n");exit(1);
      /*
      pars[1] = TruncatedMean(htmp, lowfrac, highfrac, &(pars[2]), &(errs[1]));
      errs[2]=0;
      chi[0] = chi[1] = 0;
      */
    }

    pars[2]=TMath::Abs(pars[2]);
    //hnor->SetBinContent(ix, htmp->GetBinContent(htmp->GetMaximumBin()));//htmp->Integral(0,htmp->GetNbinsX()+1));
    hmpv->SetBinContent(ix, pars[1]);
    hmpv->SetBinError(  ix, errs[1]);

    hwid->SetBinContent(ix, pars[2]);
    hwid->SetBinError(  ix, errs[2]);

    hres->SetBinContent(ix, fabs(pars[1])>EPSILON? pars[2]/fabs(pars[1]):0);
    hres->SetBinError(  ix, fabs(pars[1])>EPSILON? errs[2]/fabs(pars[1]):0);

    hchi->SetBinContent(ix, chi[1]>=1 ? chi[0]/chi[1]: 0);
    hchi->SetBinError(ix, 0);

    delete htmp;
  }

  TH1 *hhs[]={hnor, hmpv, hwid, hres, hchi};
  const TString yt[]={"N", "MPV", "#sigma", "#sigma/MPV", "#chi^{2}/NDOF"};
  const Int_t nh = sizeof(hhs)/sizeof(TH1*);
  for(Int_t ii=0; ii<nh; ii++){
    hhs[ii]->SetYTitle(Form("%s of %s", yt[ii].Data(), hh->GetYaxis()->GetTitle()));
    hhs[ii]->SetXTitle(hh->GetXaxis()->GetTitle());
    hhs[ii]->GetYaxis()->SetTitleOffset(hh->GetYaxis()->GetTitleOffset());
    hhs[ii]->SetTitle(hh->GetTitle());
  }
}

//===================================================================
//                           IO
//===================================================================

TChain * BaseUtils::InputFiles(const TString file, const TString tr, Char_t *dir)
{
  TChain *ch=new TChain(tr);
  
  if(file.Contains(".root"))
    ch->Add(file);
  else{
    ifstream fin(file);
    if(!fin){
      printf("BaseUtils::InputFiles file not found \n%s\n\n",file.Data()); exit(1);
    }
    
    TString buff;
    while(fin.good()){
      fin>>buff;
      if(buff!=""){
        if(dir){
          buff.Prepend(dir);
        }
        ch->Add(buff);
      }
    }
  }

  //const Int_t ent = ch->GetEntries(); //takes infinity time!!
  printf("\t%d trees!\n",ch->GetNtrees());
  
  return ch;
}

#endif
